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Abstract 


We  consider  error  propagation  near  an  unstable  equilibrium  state  (classified  as  an  unstable 
focus)  for  spatially  uncorrelated  and  correlated  finite-amplitude  initial  perturbations  using  short- 
(up  to  several  weeks)  and  intennediate  (up  to  two  months)  range  forecast  ensembles  produced  by  a 
barotropic  regional  ocean  model.  An  ensemble  of  initial  perturbations  is  generated  by  the  Latin 
Hypercube  design  strategy,  and  its  optimal  size  is  estimated  through  the  Kullback  -  Liebler  distance 
(the  relative  entropy).  Although  the  ocean  model  is  simple,  the  prediction  error  (PE)  demonstrates 
non-trivial  behavior  similar  to  that  existing  in  3D  ocean  circulation  models.  In  particular,  in  the 
limit  of  zero  horizontal  viscosity,  the  PE  at  first  decays  with  time  for  all  scales  due  to  dissipation 
caused  by  nonlinear  bottom  friction,  and  then  grows  faster  than  [quasij-exponentially.  Statistics  of  a 
prediction  time  scale  [the  irreversible  predictability  time  (IPT)]  quickly  depart  from  Gaussian  (the 
linear  predictability  regime)  and  becomes  Weibullian  (the  non-linear  predictability  regime)  as 
amplitude  of  initial  perturbations  grows.  A  transition  from  linear  to  non-linear  predictability  is 
clearly  detected  by  the  specific  behavior  of  IPT  variance.  A  new  analytical  formula  for  the  model 
predictability  horizon  is  introduced  and  applied  to  estimate  the  limit  of  predictability  for  the  ocean 
model. 

Keywords:  oceanography,  wind-driven  circulation,  current  prediction,  stochastic  stability, 
statistical  analysis 
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1.  Introduction 


A  number  of  recent  theoretical  studies  have  demonstrated  that  the  robust  dynamical  regimes 
[attractors]  of  oceanic  circulation  often  present  a  combination  of  [quasi-stable]  equilibrium  states  and 
transient  dynamics  between  these  flow  configurations  (Berloff  and  McWilliams,  1999;  Schmeits  and 
Dijkstra,  2000;  Stanev  and  Staneva,  2000;  Sura  et  ah,  2001;  Lermusiaux  et  ah,  2006  among  others). 
For  example,  Eremeev  et  al.  (1992)  and  Stanev  and  Staneva  (2000)  identified  single,  double,  and 
multiple  basin-scale  current  gyres  observed  in  the  Black  Sea  as  quasi-stable  equilibrium  states,  for 
which  transient  dynamics  were  induced  by  baroclinic  instability  and  mesoscale  anti-cyclonic  eddy 
activity.  Another  example  that  was  well  documented  in  observations  and  numerical  models  is  the  path 
variation  of  the  Kuroshio  south  of  Japan  (Masuda  et  al.,  1999). 

Clearly,  an  equilibrium  state  influences  the  phase-spatial  organization  of  the  local  prediction  error 
growth  rate,  i.e.  it  organizes  the  local  predictability  for  a  forecast  model.  Therefore,  the  knowledge  on 
how  small-  and  finite-amplitude  perturbations  evolve  near  this  equilibrium  state  is  important  for 
understanding  regional  model  predictability  and  identifying  the  persistence  of  circulation  ocean  and 
atmospheric  patterns  with  oscillations  near  equilibrium  states  (Robinson  et  al.,  1996). 

The  primary  goals  of  the  proposed  study  are  (1)  to  understand  what  mechanism[s]  can  fonn  the 
statistics  of  finite-amplitude  prediction  error  (PE)  near  an  unstable  equilibrium  state  identified  as  an 
unstable  focus,  (2)  to  check  how  quickly  such  statistics  depart  from  Gaussian  (if  such  a  departure  exists) 
for  the  short-  (up  to  several  weeks)  and  intermediate  (up  to  a  couple  of  months)  range  forecasts,  and  (3) 
how  to  quantify  PE  statistics  for  perfect  models  with  initial  conditions  corrupted  by  finite-amplitude 
stochastic  perturbations. 

The  computations  presented  below  assume  a  perfect  model  scenario  with  stochastic  perturbations 
in  initial  conditions.  A  non-linear  baro tropic  model  of  wind-driven  circulation  in  an  idealized  basin  is 
used  to  understand  evolution  of  prediction  error.  Although  this  model  seems  to  be  too  simple  in 
comparison  with  large  state  of  the  art  oceanic  models,  it  describes  a  generic  system  with  many  degrees 
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of  freedom  while  not  requiring  large  computer  resources.  The  latter  feature  allows  us  to  generate  non- 
rank  deficient  forecast  ensembles  and  strongly  reduce  errors  in  determination  of  appropriate 
distributions  in  comparison  with  small  ensemble  integration. 

Predictability  in  real-time  systems  can  be  qualified  in  many  ways  (for  example,  see  Smith  et  ah, 
1999).  The  proposed  study  defines  model  predictability  in  the  stochastic  stability  context.  The 
stochastic  stability  addresses  effects  of  random  perturbations  on  trajectories  of  a  dynamical  system  and 
estimates  its  stability  in  terms  of  probabilistic  measures,  such  as  expected  values  or  distribution 
functions  (Freidlin  and  Wentzell  ,1998).  In  general,  the  stochastic  stability  and  predictability  differ  from 
one  another.  However,  if  a  time  scale  quantifies  the  model  predictability,  and  if  this  scale  indicates  the 
time  when  the  forecast  uncertainty  exceeds  some  boundary  or  when  information  on  the  initial  condition 
is  lost,  the  stochastic  stability  and  predictability  are  interchangeable.  Since  these  time  scales  are  widely 
used  in  meteorology  (see,  for  example,  Toth,  1991)  and  oceanography  (Robinson  et  ah,  1996),  the 
stochastic  stability  concept  seems  to  be  a  useful  tool  for  the  predictability  analysis  of  large 
hydrodynamic  models. 

The  loss  of  superposition  and  the  extreme  inhomogeneity  common  in  nonlinear  hydrodynamic 
models  require  applying  local  measures  of  predictability  and  corresponding  time  scales  (see  Lorenz, 
1965;  Benzi  and  Camavale,  1989;  Ivanov  et  ah,  1994;  Boffetta  et  ah,  1998;  Smith  et  ah,  1999;  Mu  et  al., 
2004,  and  references  thereof).  It  is  widely  held  that  time  scales  are  related  to  the  inverse  of  the  largest 
Lyapunov  exponent  estimated  by  the  tangent  linear  models  in  assumption  of  small-amplitude  initial 
perturbations.  The  linear  approach  gives  reasonable  estimations  of  model  predictability  in  many 
practical  cases.  However,  it  cannot  provide  critical  boundaries  on  finite-amplitude  stability  of  the 
thennohaline  ocean  circulation. 

In  regional  ocean  modeling,  neither  initial  perturbations  nor  prediction  errors  are  small.  Therefore, 
the  linear  predictability  regime,  where  the  PE  grows  [quasi]  exponentially  or  even  faster  than 
exponentially,  may  be  quickly  replaced  by  the  nonlinear  predictability  regime,  for  which  the 
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predictability  time  is  much  larger  than  the  inverse  of  the  leading  Lyapunov  exponent  (Aurell  et  ah,  1996; 
Lorenz,  2005  and  others). 

To  quantify  both  the  linear  and  non-linear  predictability  the  proposed  study  uses  the  so-called 
irreversible  predictability  time  (IPT)  originally  introduced  by  Ivanov  et  al.  (1994).  Chu  et  al.  (2002) 
demonstrated  the  capability  of  IPT  for  the  analytical  estimate.  They  have  obtained  an  analytical  formula 
for  the  mean  IPT  and  variance  of  IPT  for  Lorenz-84  atmosphere  model  (Lorenz,  1984).  Here,  we  would 
like  to  show  how  IPT  statistics  change  with  a  transition  from  linear  to  non-linear  predictability  regimes, 
and  to  find  what  kind  of  PE  statistics  may  accompany  the  nonlinear  predictability  regime. 

The  paper  is  organized  as  follows.  Section  2  briefly  discusses  IPT  and  its  statistics.  Section  3 
analyzes  the  specificities  of  the  reference  solution  reproduced  by  a  barotropic  regional  ocean  model. 
The  ensemble  of  stochastic  perturbations  added  to  the  initial  conditions  is  described  in  Section  4. 
Section  5  depicts  the  optimal  ensemble  size  of  stochastic  perturbations  using  the  Kullback-Leibler 
distance.  Section  6  discusses  the  initial  PE  decay  due  to  viscosity  damping.  Section  7  studies  the 
response  of  the  model  to  finite-amplitude  initial  perturbations,  and  analyzes  a  basic  feature  of  the  PDF 
of  IPT  (denoted  r-PDF),  such  as  non-Gaussianity.  Section  8  gives  evidence  for  using  a  three- 
parameter  Weibull  distribution  as  IPT  statistics  in  the  non-linear  predictability  regime.  Variations  of 
mean  IPT  and  its  variance  for  transition  from  the  linear  to  nonlinear  predictability  regime  are  discussed 
in  Section  9.  Section  10  estimates  the  model  predictability  horizon  (the  maximum  predictability  time  for 
the  given  statistics  of  initial  perturbations).  Section  1 1  summarizes  the  obtained  results. 

2.  Irreversible  Predictability  Time 

Let  the  prediction  error  Z(x,  t)  be  defined  as  a  difference  between  the  reference  solution  Y(x,t) 
[i.e.,  the  solution  of  a  perfect  model  without  errors  in  initial  conditions  (Lacarra  and  Talagrand,  1988)] 
and  an  individual  forecast  Y(x,t) : 

Z(x,t)  =  Y(x,t)-Y(x,0,  Z0  =Z  (x,f0), 


5 


where  (x,  t)  are  spatial  coordinates  and  time;  Z0  refers  to  the  initial  perturbations.  Y (x,t)  is  the  state 

vector,  which  may  include  velocity,  temperature,  salinity  and  other  fields. 

The  ocean-atmospheric  model  predictability  is  often  quantified  by  the  weighted  relative  root  mean 
square  error  (see  Robinson  and  Haidvogel,  1980;  Holland  and  Malanotte-Rizzoli,  1989;  Brasseur  et  ah, 
1996,  Robinson  et  ah,  1996;  Wirth  and  Ghil,  2000  among  others)  written  as 

J(Z0 ,  w,  0  =  (z, wz)  •  J~lorm  (t),  ( 1 ) 

where  W  is  the  weight  matrix,  is  the  inner  product ,  the  function  Jnorm  (/)  is  specified  from  the 

physics. 

The  IPT  is  defined  as  the  time  r  ,  at  which  J(Zq,W  ,t)  reaches  a  predetermined  level  s~  for  the 
first  time: 

r(z0,W,j)  =  inf(f  j(z0, W,t)>  £21,  (2) 

t>  ov  J 

where  s  is  a  non-dimensional  tolerance  level  (accepted  prediction  accuracy). 

This  definition  is  illustrated  by  Fig.  la.  Clearly,  the  IPT  defines  the  model  predictability  on  the 
condition  that  any  returns  of  model  predictability  (the  shaded  zones  in  Fig.  la)  do  not  contribute  to  the 
prediction  skill. 

The  mean  IPT  differs  from  the  e-folding  or  the  doubling  time  when  J  oscillates  or  is  random.  To 
compare,  for  example,  the  e-folding  and  the  irreversible  predictability  time,  we  should  suppose  in  Eq. 

(2)  Jnorm  =  (Zq,Zq)  and  s2  -  e1 .  The  e-folding  time  is  the  time  when(/)  crosses  e 2  (Fig.  lb), 

re  (W)  =  maxjj  |  (j( Z0 ,  W, t))  <  e2  J ,  (3) 

where  the  brackets  (...)  denote  the  average  over  the  ensemble  of  initial  perturbations  Z0  . 

The  mean  IPT  for  the  same  e  is  computed  by 

(r(W,e))  =  ^inf(f  |  j(z0,  W,t)  >  e2Jj  ,  (4) 
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where  the  averaging  is  over  the  ensemble  of  IPTs  (rl,...,rN)  induced  by  the  ensemble  of  Z0  (Fig.  1  c). 
Chu  and  Ivanov  (2005)  pointed  out  that  the  mean  IPT  is  the  lower  bound  of  e-folding  time: 

re(W)  >  (r(W,e)). 

In  practical  applications,  the  ensemble  generated  z  -  PDF  [  F(z,  W,  s) ,  hereafter,  this  is  a  PDF 
which  corresponds  to  the  given  ensemble  of  initial  perturbations],  the  cumulative  distribution  function  of 
r  (r-  CDF) 

P(W,s,t-t0)  =  Probfy  >  t-t0)  ,  (5a) 

and  r  -  moments  calculated  by 

00 

(zk(W,£))  =  kj(t-t0)k~lP(W,£,t-t0)dt,  k=l...K,  (5b) 

A 

may  be  used  to  quantify  model  predictability.  The  first  four  unbiased  r  -  moments  detennine  z  -  mean 
((r)),  r-  variance  ((^Sz2^),  r-  skewness  (SK)  and  z-  kurtosis  (KU). 

For  simplicity,  the  further  analysis  supposes  (a)  to  replace  the  weight  matrix  W  by  the  identity 
matrix  I ,  (b)  to  describe  flow  dynamics  in  a  quasi-geostrophic  approximation  using  the  geostrophic 
stream  function  P  (Pedlosky  1987) ,  and  (c)  to  take  J norm  =  (P,P) . 

3.  The  Reference  Solution 

A  shallow  water  circulation  computed  in  a  flat  bottom  semi-enclosed  basin  and  forced  by  wind 
and  water  flux  across  its  open  boundary,  is  taken  as  the  reference  solution.  Our  model  leaves  out  the 
effects  of  topographic  and  baroclinic  processes  but  it  reproduces  a  highly  non-linear  flow  with  a  balance 
on  the  ft  -plane  between  steady  wind  forcing,  nonlinear  bottom  friction  and  inertial  terms  of  the  model. 
The  model  is  similar  to  Veronis’  model  of  wind-driven  circulation  in  a  rectangular  basin  (Veronis, 
1966)  but  with  a  non-linear  bottom  friction. 
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The  computation  domain  presents  a  rectangular  basin  with  depth  H  (=  2  km)  centered  at  35°N. 
The  domain  is  bounded  by  the  rigid  ( A )  and  open  boundaries  (Domain-A  depicted  in  the  top  panel  of 
Fig.  2).  This  basin  extends  L2  =  1000  km  [Lx  =  1050  km]  in  the  north-south  (east-west)  directions.  A 
Cartesian  coordinate  system  is  used  with  the  origin  in  the  southwest  comer.  The  x1  -axis  points  towards 
the  east,  and  the  x2  -axis  towards  the  north. 

The  barotropic  mode  of  the  Princeton  Ocean  Model  (POM)  (Blumberg  and  Mellor,  1987)  is 
applied  to  compute  the  horizontal  velocities  ux  ,u2 ,  and  the  surface  elevation  /  in  Domain-A  with  no¬ 
slip  boundary  conditions.  The  circulation  is  forced  by  wind  with  stress  varying  with  latitude 

—  =  (~10~3w2 /5|2)cos  7TX'1  ,  (6a) 

A)  l2 

•5 

and  the  prescribed  open  boundary  conditions  (ub and<^)  explained  below.  Here,  p0  =  1025.44  kg/nr 

is  the  reference  density.  The  model  runs  with  time  step  equaled  to  2.5  minutes  and  reproduces  the 
circulation  with  horizontal  resolution  of  50  km. 

The  Coriolis  parameter  varies  linearly  with  a  beta  plane  approximation  /  =  f0  +  J3  x2,  where 
/o  =  2/2 sin  cp()  and  ft  =  (2/2/  a)  cos  cp0  .  Here,  Q  and  a  are  the  rate  of  rotation  and  the  radius  of  the 

Earth,  respectively;  cpQ  =35°;  f0  =  0.73  x  10~4  s"1 ,/?  =  2x  10_1 1  m'1^-1 . 

The  horizontal  kinematic  viscosity  is  set  to  zero.  The  bottom  stress  is  parameterized  by  the 
quadratic  drag  relation: 


where  the  drag  coefficient  a  =0.0025.  No  model  spin  up  exists  for  a  <0.0025  when  any  solution  is 
unstable. 

The  open  boundary  conditions  are  specified  by  Chu  et  al.’s  (1997)  approach.  Accordingly  to  this 
approach  POM  is  firstly  integrated  with  wind  stress  (6a)  and  dissipation  (6b)  from  rest  ( u  =  v  =0)  and 
flat  surface  (^  =0)  for  150  days  in  a  closed  rectangular  basin  fonned  by  extension  of  Domain-A  up  to 
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2000  km  in  the  east  (Domain-B  shown  in  the  bottom  panel  of  Fig.  2).  The  velocity  and  surface  elevation 
in  the  left  part  of  Domain-B  computed  on  day- 10  are  utilized  as  the  initial  conditions  for  the  control  run 
in  Domain-A.  The  values  of  ux  and  C,  along  the  line  x,  =  1 050  km  are  the  prescribed  open  boundary 
conditions  ub  and  C,h  between  day-10  and  day-140.  The  model  spin-up  was  found  in  Domain-A  after  60 

day  integration.  Here,  the  kinetic  energy  of  the  reference  flow  slowly  decays  with  the  power  exponent 
of  1000  day'1,  oscillating  with  period  of  about  180  days. 

The  circulation  pattern  evolves  from  a  single  semi-closed  gyre  with  a  maximum  velocity  of 
about  0.35  m  s'1  (Fig.  3a)  and  sea  surface  elevation  between  0.05  to  0.1  m  (not  shown)  to  a  multi-gyre 
structure  with  maximum  velocities  up  to  0.9- 1.0  m  s'1  (Fig.  3b)  and  high  surface  elevation  near  1  m  in 
the  west-northern  part  of  the  basin  (not  shown).  This  multi-gyre  structure  is  a  model  spin-up  reached 
after  60  day  integration. 

The  model  spin  up  is  a  nonlinear  regime  identified  by  two  non-dimensional  numbers  (Pedlosky, 
1987): 

Rc>\  =  max(-^— )  ~  0.7  -0.1,  Ro2  =  max(— *— )  ~  0.12,  (7) 

foL  f0T 

where  f0~  7  x  10'5  s'1  is  the  Coriolis  parameter;  U~  0.5-0. 7  m  s'1  is  the  characteristic  velocity  in  the 

basin;  L  ~  105  m  and  T  ~  1  day  are  the  characteristic  spatial  and  temporal  scales  of  the  flow, 
respectively. 

The  chosen  model  configuration  allows  us  to  use  the  same  model  for  the  analysis  of  model 
predictability  affected  by  different  sources  of  uncertainty:  stochastic  errors  inserted  in  initial  condition 
(the  present  study),  wind  (Ivanov  and  Chu,  2006),  and  open  boundary  conditions  (Ivanov  and  Chu, 
2007).  Cross-correlations  between  different  types  of  errors  can  also  be  studied. 

Using  the  quasi-geostrophic  approximation  allows  introducing  a  model  phase  space  with  the  basis 
composed  from  M  eigen-functions  [ y/i,...,y/ M  ]  of  a  plane  Laplacian  operator  (Aj_ )  defined  in  Domain- 
A  (for  details  see  Appendix  A)  and  get  the  following  spectral  representations  for  the  reference  solution 
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(8a) 


M 

*)P(x,  t)  —  Am  {t)y/  m  (x)  +  (x,  t ) , 

m=l 

and  a  perturbed  solution 

M 

S'n*,t)  =  Y SAm  (OWm  (x) >  (8b) 

m=\ 

where  SlF  =  W  -  '/y ,  IP  is  an  individual  forecast,  *P/t(2rm  is  the  hannonic  function  described  in 
Appendix  A. 

The  spectral  coefficients  [Afyt),...,  AM(t)]  and  [A4i(t),...,A4^(t)]  define  the  reference  and 
perturbed  trajectories,  starting,  respectively,  from  points  [A® ]  and  [SA® SA^  ]  in  the  phase 
space.  Appendix  A  specifies  the  optimal  truncation  number  for  both  the  reference  and  perturbed  flows 
as  M  «  102 . 

Therefore,  the  dimensionality  of  a  dynamical  system  ( M  )  embedding  the  ocean  model  should  be 
-0(10  ).  However,  the  real  number  of  degrees  of  freedom  (M0)  determined  by  Syrovich’s  approach 
(Syrovich,  1989;  Aubry  et  al.  ,1991)  is  only  about  15  (Chu  and  Ivanov,  2005). 

We  used  a  well-known  classical  stability  analysis  (Guckenheimer  and  Holmes,  1983)  to  classify 
the  model  spin  up  (quasi-equilibrium  state)  in  dynamical  system  context.  The  governing  equations  in  the 
spectral  form  (not  shown)  were  linearized  near  the  spin  up  state,  and  eigenvalues  and  eigenfunctions  of 
appropriate  tangent  operator  were  calculated.  The  spectrum  of  the  eigenvalues  contained  positive  and 
negative  real  numbers  only.  Thus,  the  quasi-equilibrium  state  is  an  unstable  focus  in  the  model  phase 
space.  Any  model  trajectory  tends  to  move  toward  (off  from)  this  point  in  a  phase  plane  determined  by  a 
pair  of  low-order  (high-order)  modes. 

4.  Perturbations  of  Initial  Conditions 

The  main  question  in  ensemble  forecasting  is  how  to  generate  a  set  of  initial  perturbations.  The 
present  study  uses  a  Monte-Carlo  method  to  produce  ensembles  of  perturbed  solutions  from  running 
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multiple  POM  simulations  with  initial  perturbations,  which  are  drawn  randomly  from  a  specified 
probability  density  function. 

It  is  often  assumed  in  ensemble  forecasting  that  for  a  forecast  model  of  nearly  I  05  degrees  of 
freedom,  direct  sampling  of  small  size  [~0(1O  )]  ensembles  is  of  limited  utility  because  the  inherent 
sampling  error  may  overwhelm  the  desired  covariance  infonnation.  To  correct  these  distortions  different 
selective  sampling  procedures  are  used.  See  Palmer  (2000)  for  extensive  review  of  different  techniques 
of  ensemble  generation  including  those  based  on  singular  (Moltenti  and  Palmer,  1993)  and  breeding 
vectors  (Toth  and  Kalnay,  1997).  Miller  and  Enret  (2002)  examined  small  ensembles  drawn  from 
spaces  spanned  by  singular  vectors  and  by  bred  vectors  for  nonlinear  dynamical  systems  with  multiple 
attractors. 

A  sound  approach  for  estimating  model  predictability  was  developed  by  Lermusiaux  and  Robinson 
(1999),  Lennusiaux  et  al.  (2000),  and  Lermusiaux  (2001,2002).  The  approach  called  the  “error  subspace 
statistical  estimation”(Lermusiaux  and  Robinson,  1999)  combines  the  dynamical  equations  for  the  ocean 
state  in  their  discrete  form  in  space,  with  a  decomposition  of  error  covariance  to  initialize  and  evolve  the 
“dominant”  eigendecomposition  of  the  variability  covariance  matrix,  merging  data  and  dynamics.  An 
essential  feature  of  the  dominant  eigenvectors  is  that  they  indicate,  evolve  and  organize  the  directions  in 
the  variable  space  that  have  largest  statistical  significance,  based  on  a  variance  measure.  Therefore,  the 
approach  provides  a  framework  for  investigation  of  large  and  complex  system  like  the  ocean. 
Lermusiaux  et  al.  (2006)  have  demonstrated  explicit  capability  of  the  approach  through  a  number  of 
practical  examples  including  acoustical  and  biological  models. 

Lor  a  number  of  reasons  we  simply  sample  initial  perturbations  from  specified,  multivariate 
nonnal  distributions.  Lirst,  due  to  simplicity,  our  hydrodynamic  model  is  able  to  produce  forecast 

ensembles  containing  up  to  ~O(105)  perturbations.  Such  large  ensembles  have  desirable  statistical 
characteristics  as  revealed  by  a  series  of  specialized  statistical  tests  considered  in  Section  5. 
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Second,  the  Latin  hypercube  (LHC)  design  (Latin  Hypercube,  2001)  was  applied  to  simulate  a 
highly  unifonn  distribution  of  an  initial  error  in  the  phase  space.  Using  pure  probabilistic  arguments, 
Downing  et  al.  (1985)  pointed  out  that  the  Latin  hypercube  design  is  more  effective  than  the  classical 
Monte-Carlo  method.  For  obtaining  dense  error  coverage  with  the  same  degree  of  homogeneity  from 

Monte-Carlo  samples  and  through  the  LHC  design  NM°  and  N  ■  (2Mq  +  2)  statistical  realizations  are 
required,  respectively.  Here,  N  is  the  number  of  statistical  realizations  necessary  to  simulate  one  degree 
of  freedom.  Typically,  N  ~  0(1 03 ) . 

For  Mq  =15,  the  classical  Monte-Carlo  method  should  require  ~O(10 45  )  initial  perturbations 
for  a  statistically  significant  estimate.  This  is  not  feasible  with  available  computer  resources. 
Comparable  results  can  be  obtained  by  the  LHC  design  approach  with  only  ~O(104)  initial 
perturbations. 

Third,  our  study  focuses  only  on  the  physics  of  finite-amplitude  errors  and  their  contributions  to 
losing  model  predictability.  Thus,  small  forecast  ensembles  similar  to  those  utilized  by  operational 
forecast  models,  are  not  considered  here. 

The  initial  perturbations  u'  =  [«[  (x),«2  (x)  ]  are  assumed  to  be  2D  isotropic  Gaussian  white  noise 
[white  noise-like  perturbations  (WNLPs)]  with  the  two-point  correlation  function 

(u'i(x)u'j(x'))  =  I2SijS(x-x'),  (9) 

or  2D  isotropic  Gaussian  spatially  correlated  noise  [red  noise-like  perturbations  (RNLPs)]  with  the 
two-point  correlation  function 

(u'i ( x)u'j (x'))  =  I2Sj  exp  -  **  ,  (10) 

2 

where  x  =  (jc15jc2),  ( R± ,  /  )  are  correlation  radius  and  noise  variance  (intensity  of  perturbations), 


respectively.  These  perturbations  are  directly  added  to  the  initial  conditions.  The  technical  details  of 
generating  Gaussian  noises  with  correlation  functions  (9)  and  (10)  can  be  found  in  Sabel’feld  (1991). 
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—2  2  2  2  2  —2 

The  non-dimensional  intensity  of  the  initial  perturbations  /  =/  / 7q  (7q  =  I  irT.v  )  will  be  used 
below. 

Noise  with  correlation  functions  (9)  and  (10)  is  the  popular  model  of  errors  for  optimal 
interpolation  or  spline  fitting.  Both  of  these  procedures  are  applied  to  construct  initial  conditions  for 
ocean  models  from  irregularly  spaced  data  (Brasseur  et  ah,  1996;  Robinson  et  ah,  1996  and  others). 

5.  Optimal  Ensemble  Size 

The  LHC  design  approach  provides  the  dense  error  coverage  of  the  model  phase  space  for 
~O(104)  initial  perturbations.  This  number  of  initial  perturbations  is  a  trade-off  between  the  ensemble 
ability  to  reproduce  main  features  of  PE  statistics,  and  the  computational  cost.  However,  the  optimal 
number  of  initial  perturbations  (Nopt)  should  be  specified  for  the  concrete  ocean  model.  We  suggest  to 

estimate  this  number  through  the  Kullback-Leibler  (KL)  distance  (the  relative  entropy)  (White,  1994). 

The  KL  distance  is  a  natural  distance  function  from  a  ’’true”  probability  density,  F^ ,  to  a  “target” 
probability  density,  FN .  For  continuous  density  functions,  the  KL  distance  is  defined  as 


uu 

KL  N  (fx  \Fn  )  =  J  d  t  Fn  (r)  log 


^oo  (r) 
Fn(t)' 


(ID 


where  Fn(t)  and  F-fj{z)  are  z  -  PDFs  computed  for  an  N  sample  ensemble  and  a  hypothetic  ensemble 
with  infinite  sampling,  respectively.  In  practice,  a  difference  between  two  distributions  is  negligible  if 
KLn  <  5.0x1  0“3  (White,  1994). 

To  calculate  the  KL  distance  we  suppose  Fx  =  F\  ooooo  because  only  small  differences  in 

r  -  statistics  estimated  from  ensembles  of  5.0  x  103 ,1.0x  104  ,  2.0  x  104 ,5.0x  104  and  1.0  x  105  samples 
have  been  observed.  Typical  behavior  of  KL ^  with  the  growth  of  N  is  shown  in  Fig.  4a.  The  KL ^ 
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rapidly  reduces  with  N  from  2.0x10  1  (IV  =20)  to  4.0x10  3  ( N  =  103).  Therefore,  Nopt  was  chosen 
as  103 . 

3 

The  first  four  moments  of  IPT  also  quickly  converge  as  N  increases  to  10  :  r-mean  to  43.8  day 
(Fig.4b),  r  -  most  probable  to  42.6  day  (Fig.4c),r -variance  to  8.8  day  (Fig.  4d),  and  z  -  skewness  to 
0.75  (Fig.  4e).  Kurtosis  is  most  sensitive  to  the  choice  of  N  (Fig.  4g),  and  varies  between  3.9  and  4.1 
for  N~10  .  However,  these  variations  are  quite  small  to  mask  the  non-Gaussian  feature  of  IPT  statistics 
(Gaussian  statistics  is  identified  by  both  SK=0.0  and  KU=  3.0). 

By  varying  characteristics  of  initial  perturbations  (/  and  R]_)  we  have  found  that  Nopt  =10  is  an 

acceptable  choice  for  any  combinations  of  /  and  R]_ .  Thus,  all  forecast  ensembles  below  have  large 
size  ( N opt  =  1 0  )  and  are  not  rank-deficient  (Nopt»M0).  This  reduces  the  sampling  errors 
significantly. 

6.  Prediction  Error  Evolution 

6.1  Different  stages  of  PE  evolution 

A  number  of  stages  for  PE  evolution  (at  least  four)  are  observed  for  both  WNLPs  (Fig.  5a)  and 

RNLPs  (Fig.  5  b).  All  these  stages  are  clearly  identified  by  the  growth  rate  Q  =  —  In  <  J  >  (Fig.  6a, b). 

dt 

Initial  error  decay  is  observed  for  the  first  ten  days  of  PE  evolution  (Fig.  5a, b)  where  the  growth 
rate  Q  evolves  as 

0~0o[exP(«oO-«i],  (12) 

where  is  the  decay  exponent,  Q§a\  ~  -0.45  ,  ag  =  In  a\lt\,  t\  =10  days. 

Non-exponential  initial  error  decay  corresponding  to  (12)  differs  from  a  quasi-exponential  decay 
obtained  in  Wirth  and  Ghil’s  (2000)  model  with  the  dissipative  operator  v  Aj_  and  quite  large  horizontal 


14 


viscosity  v  .  In  contrast  to  this  model  the  zero  horizontal  viscosity  is  used  in  our  model.  Therefore,  the 
non-exponential  decay  seems  to  be  caused  by  nonlinear  bottom  friction. 

Then  (J(t))  remains  quasi-stationary  low  values  during  10  and  20  days  for  WNLPs  (Fig.  5a) 
and  RNLPs  (Fig.  5b),  respectively.  At  this  stage  (stagnation  stage)  the  growth  rate  Q  oscillates  near 
zero  (Fig.  6a, b). 

During  the  third  stage  (after  day-20  and  day-30  for  WNLPs  and  RNLPs,  respectively)  PE  grows 
faster  than  exponentially  (Fig.  5a, b).  Fig.  6a, b  shows  that  here  the  growth  rate  increases  with  t  linearly, 

Q~t  (13) 

Mechanisms  for  PE  to  grow  faster  than  [quasij-exponentially,  were  discussed  earlier  in  the 
scientific  literature  (Lacarra  and  Talagrand,  1988;  Smith  et  ah,  1999  among  others).  Therefore,  they  are 
not  in  the  focus  of  the  present  study.  Duration  of  “super-exponential”  error  growth  does  not  exceed  7-10 
days. 

Nonlinear  interactions  among  various  scales  slow  down  (Fig.  6  a,b),  and  then,  limit  further  growth 
of  PE  (Fig.  5a, b).  The  transition  from  the  linear  to  nonlinear  predictability  regime  is  identified  by  a 
behavior  of  the  growth  ratio:  Q  reaches  a  maximum  value  and  then  quickly  reduces.  Although  linear 
effects  can  also  stimulate  similar  behavior  of  Q  (Smith  et  ah,  1999),  more  detailed  analysis  provided  in 
Section  7  gives  evidence  for  the  nonlinear  predictability  regime. 

The  duration  of  the  nonlinear  predictability  regime  (up  to  10-15  days)  does  not  seem  to  be  short, 
and  is  comparable  with  the  time  when  the  PE  grew  faster  than  exponentially.  The  larger  the  error 
amplitude,  the  faster  the  transition  from  linear  to  nonlinear  predictability  regime.  Therefore,  the 
nonlinear  regime  should  play  an  important  role  in  understanding  model  predictability  and  its  limits. 

6.2  Initial  error  decay 

Initial  error  decay  was  documented  earlier  by  Lorenz  (1996),  Brasseur  et  al.  (1996),  Molteni  et 
al  (1996),  Wirth  and  Ghil  (2000),  Vannitsem  and  Nicolis  (1997),  Snyder  et  al.  (2003)  and  others.  Their 


15 


results  show  (a)  the  relation  of  the  decay  to  model  dissipation,  and  (b)  that  stability  of  the  reference 
solution  estimated  by  the  kinetic  energy  norm  is  less  than  when  other  norms  are  used. 

The  initial  error  decay  quantified  by  the  kinetic  energy  norm  diminishes  or  even  disappears  with 
reducing  the  coefficient  of  horizontal  viscosity  v/7 ,  when  dissipation  process  were  parameterized  by  a 
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hyper-dissipation  operator  -n/7  Aj_  (Snyder  et  al.,  2003;  Vannitsem  and  Nicolis,  1997). 

In  contrast  to  this,  Wirth  and  Ghil  (2000)  have  demonstrated  that  the  square  root  of  the  stream- 
function  error  variance  ( J )  for  randomly  inserted  small-amplitude  perturbations  first  decays  due  to 
viscous  damping  parameterized  by  the  usual  dissipative  operator  vA±.  Although  for  smaller  values  of 
horizontal  viscosity  v  the  decay  was  slower,  it  did  not  disappear  with  reducing  v  . 

Our  computations  demonstrate  the  existence  of  the  initial  error  decay  for  both  the  kinetic  energy 
nonn  and  norm  (1)  only  due  to  nonlinear  bottom  friction  if  the  friction  coefficient  a  >  a  .  The  choice  of 
J  norm = '  'n  norm  (1)  does  not  change  the  obtained  result,  although  the  decay  becomes  slower.  The 
decay  also  exists  for  non-zero  coefficients  (v  and  i//7)  of  horizontal  viscosity. 

Since  the  initial  error  decay,  in  general,  does  not  depend  on  the  amplitude  of  initial  perturbations 
(Fig.  5  a,  b),  to  understand  its  physics  we  limit  ourselves  to  the  case  of  small-amplitude  initial 
perturbations.  In  our  opinion,  the  initial  error  decay  may  be  explained  by  features  of  a  nonlinear 
dissipation  scheme  and  statistical  features  of  initial  perturbations  with  correlation  functions  (9)  and  (10). 

Using  the  assumption  on  small  amplitudes  of  initial  perturbations,  the  term,  which  is  responsible  for 
PE  dissipation  (see  Appendix  B),  can  be  written  as: 


F\u\ ) «  aE 


-1/2 


/l\(u'l)2)  +  /2  {ulu2  ) 


F2u'2)*aE~V2 


rAu'2?)  +  r2(«i'«2). 


(14) 

(15) 
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where  ( ul,u2 )  and  ( u[,u'2 )  are  velocities  of  the  reference  and  perturbed  solutions,  respectively; 
E  =  w i 2  +  w 2 2  is  the  kinetic  energy  of  the  reference  flow,  y{  =  2  /7,2  +  w22 ,  y2=ulu2  and 
7s  =  «i2  +2w22. 

During  the  first  stage  (initial  error  decay),  the  PE  dissipation  is  only  determined  by  the  first  terms 
in  the  right  hand  sides  of  Eqs  (14)  and  (15),  because  the  initial  perturbations  are  statistically  isotropic: 

(u,i(x)u'j(x)\  =  0 ,  if  i  *  j  (the  isotropy  feature).  Both  renormalized  coefficients  of  friction  a  E~^  2 iq 

and  a  E  yj  are  positive  and  grow  with  increasing  kinetic  energy  as  E 

During  the  stagnation  stage  (Fig.  5  a,  b),  the  correlations  between  u[  and  u2  become  not 
small.  The  growth  of  the  renormalized  coefficients  is  compensated  by  the  tenn  y2 (it [u2)<  0 .  This 
process  stops  the  PE  decay.  After  day-20  for  WNLPs  and  day-30  for  RNLPs,  y2(u[u2)  <  0  continues 
to  reduce  the  dissipative  terms  (Fxu[)  and  ( F2u'2 ),  what  reduces  the  model  dissipation  and  stimulates 
the  growth  of  prediction  error. 

Nonlinear  bottom  friction  (6b)  leads  to  effective  decay  of  the  prediction  error  at  all  spatial  scales 

larger  than  a  scale  detennined  by  the  100-th  mode  (see  Fig.  7a  and  Fig.  7b).  At  this  stage  of  PE 

evolution  geostrophic  adjustment  inducing  an  upscale  flux  of  the  prediction  error  does  not  play  an 

important  role.  A  measure  of  re-distribution  of  the  kinetic  energy  among  modes  is  the  spectral  entropy 

(Aubry  et  al.  ,1991).  The  entropy  is  equaled  to  1  when  the  kinetic  energy  is  homogeneously  distributed 

among  all  the  modes  (maximum  disorder),  and  to  0  when  the  energy  is  contained  in  a  single  mode 

(maximum  order).  The  spectral  entropy  is  computed  by 

M  M 

S  —  — (log TT)  Pm  lo§/?7?j5  Pfj]  ~bm  lb  ,  b  —  b]ri  ,  (16) 

m=m0  m=m0 

where  bm  are  Am  ( 8Am )  for  the  reference  (the  perturbed)  solution;  m0  and  M  determine  the  spectral 
band  for  calculating  the  spectral  entropy. 
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The  spectral  entropy  computed  for  the  PE  (the  reference  solution)  is  shown  in  Fig.  8a  (Fig.  8b).  Fig. 
8a  does  not  show  any  significant  re-distribution  of  kinetic  energy  among  different  flow  scales  up  to  10- 
12  days.  A  strong  upscale  flux  of  PE  is  observed  for  low-order  modes  (from  1  to  20)  after  day-12:  the 
spectral  entropy  approximately  reduces  into  half.  Such  a  process  seems  to  be  considerably  weaker  for 
high-order  modes.  In  contrast  to  the  PE,  the  complexity  of  the  reference  flow  monotonically  grows  up  to 
day-40  (Fig.  8  b). 

7.  Response  To  Finite  -Amplitude  Initial  Perturbations 


Finear  intuition  suggests  that  the  larger  the  amplitude  of  the  initial  perturbations,  the  higher  the 
probability  of  obtaining  low  model  predictability.  Prediction  errors  should  steadily  increase  with  a 
prediction  time  scale  in  the  linear  predictability  regime.  In  contrast  to  this,  forecasts  skill  may  decay 
slower  when  amplitude  of  initial  perturbations  grows.  Our  computations  show  that  the  growing 
perturbations  rapidly  adopt  a  horizontal  scale  comparable  to  that  of  the  reference  state  (linear 
predictability  regime),  and  their  further  growth  is  limited  by  interactions  with  this  state  and  among  them 
(nonlinear  predictability  regime).  In  the  nonlinear  predictability  regime  the  PE  demonstrates  clear 
contributions  from  the  cumulative  effects  of  flow  scales,  and  the  predictability  time  is  no  longer 
measured  by  the  inverse  of  the  leading  Fyapunov  exponent.  Moreover,  model  predictability  enhances 

with  the  growth  of  the  correlation  radius  R±  and  is  less  sensitive  to  the  choice  of  the  intensity  /  . 

To  understand  correlations  between  model  predictability,  the  amplitude  of  initial  perturbations,  the 
correlation  radius  and  the  noise  intensity,  let  us  show  how  the  amplitude  of  initial  perturbations  is 


affected  by  Rj_  and  /  .  Using  the  spectral  representation  for  noise  covariance  matrix  obtained  in 
Appendix  C,  the  maximum  amplitude  of  initial  perturbations  is  estimated  as 


((-SO2)''2  =(2x)'/1LiL2R11I  (17) 


where  /Lmax  is  the  maximum  eigenvalue  of  a  matrix  detennined  in  Appendix  C. 
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Eq.  (17)  clearly  shows  that  the  amplitude  of  initial  perturbations  is  most  sensitive  to  the  correlation 
radius  but  not  the  noise  intensity,  because  the  amplitude  may  grow  as  R]_  and  only  linearly  with  /  . 

Our  computations  demonstrate  the  existence  of  the  linear  predictability  regime  identified  by  non- 
Gaussian  statistics  and  [quasi]-exponential  growth  of  PE,  for  a  small  correlation  radius  (i?j_  «  50  km, 

-9 

the  case  of  WNLPs)  and  large  noise  intensities  (/  <  0.2  ).  For  WNLPs,  typical  z  -  PDF  was  close  to 

Gaussian  if  Z  —  0.01-0.05.  The  growth  of  /  up  to  0. 1-0.2  resulted  in  a  weak  asymmetry  for  the 
r-PDF  (SK— >0.15),  and  departs  from  non-Gaussian  (KU  —>3.10).  However,  although  such  a 
t  -  PDF  has  a  short  tail  (labeled  by  1  in  Fig.  9a),  it  was  still  close  to  Gaussian  and  the  mean  z  -  IPT 
reduced  with  the  growth  of  amplitude  of  initial  perturbations.  The  z  -  PDF  quickly  departs  from 

Gaussian  with  the  growth  of  /  after  0.2.  However,  from  the  physical  point  of  view,  such  initial 
perturbations  seem  to  be  too  large  to  exist  in  reality. 

The  nonlinear  predictability  regime  appears  as  R j_  grows.  Both  highly  non-Gaussian  z  -  PDFs, 
and  the  mean  IPT  that  grows  with  i?j_ ,  indicate  that  the  PE  becomes  nonlinear.  A  typical  non-Gaussian 
z  -  PDF  (SK«0.8,  KU«4.0)  computed  for  a  finite  correlation  radius  is  demonstrated  in  Fig.  9  b.  The 
long  PDF  tail  (labeled  by  2  in  Fig.  9b)  is  clearly  seen  in  this  figure.  The  tail  is  formed  by  rare  individual 
forecasts  (IPT  up  to  60  days),  each  of  which  is  longer  than  the  mean  ensemble  forecasting  (IPT  of  about 
44  days). 

Asymmetry  of  z  —  PDFs  becomes  higher  for  the  larger  values  of  correlation  radius  R  (Fig.  10 
c).  SK,  which  is  a  measure  of  asymmetry,  increases  up  to  0.8  when  R  tends  to  100  km.  Farger  values 
of  mean  IPT  (Fig.  10a)  and  z  -  variances  (Fig.  10b)  correspond  to  more  asymmetric  PDFs.  Highly  non- 
Gaussian  (KU  «  4,  Fig.  lOd)  and  sharp  z-  PDFs  with  long  tails  stretching  to  large  prediction  times 
accompany  this  nonlinear  predictability  regime.  The  explicit  growth  of  mean  predictability  time 
observed  with  the  growth  of  correlation  radius  R±  is  a  strong  evidence  of  the  nonlinear  predictability 
regimes  caused  by  inhomogeneous  morphology  of  the  model  phase  space  (Kaneko,  1998). 
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8.  Weibull  Statistics 


The  results  above  demonstrated  a  fast  departure  of  v  -  PDFs  from  Gaussian  shape  with  the  growth 
of  amplitude  of  initial  perturbations.  For  finite-amplitude  initial  perturbations  r  -  PDFs  were  highly 
non-Gaussian  and  asymmetric.  The  following  question  arises:  what  distribution  is  the  best  lit  for 
such  r- PDFs?  If  the  analytical  form  of  such  a  distribution  can  be  found,  it  can  be  useful  for  the 
parametric  estimate  of  ensemble  generated  PDFs  from  limited  observation  samples  and  small  forecast 
ensembles. 

Our  computations  demonstrate  that  the  tailed  r  -  PDFs  reconstructed  by  a  non-parametrical 
technique  based  on  the  Epanichenikov’s  kernel  and  the  bootstrap  re-sampling  procedure  (Good,  1996) 
directly  from  the  Monte-Carlo  samples,  are  fitted  by  the  three-parameter  Weibull  distribution  function 


(Bury,  1999) 


'  T-y'\ 

p-\ 

f  \ 

r  -  / 

exp 

— 

\  V  J 

l  V  J 

(18) 


Here,  ?7, /and  ft  are  scale,  shape  and  location  parameters.  We  found  no  difference  between  the 
reconstructed  PDFs  and  their  Weibull  counterparts  with  at  least  a  95%  confidence  level. 

We  suggest  to  identify  the  distribution  parameters  of  (18)  by  the  probability  weighted  moments 
(PWMs).  Definition  of  PWMs  and  their  features  are  given  in  Greenwood  et  al.  (1979).  For  a  number  of 
reasons  these  moments  seems  to  be  more  attractive  for  the  practical  calculations  of  prediction  scales 
than  the  classical  statistical  moments  (CSMs).  First,  the  PWMs  are  less  sensitive  to  sampling  than  the 
CSMs.  Second,  the  classical  statistical  moments  may  not  exist  for  PDFs  with  long  (“heavy”)  tails.  For 

example,  if  a  cumulative  distribution  function  P  has  asymptotic  tail  ~  t~a  as  t  — >  go  ,  then  the  classical 
t  -  moments  of  (k-cr-1)  order  do  not  exist  because  integral  (5)  does  not  converge.  In  contrast  to  the 
CSMs,  the  PWMs  exist  for  any  tailed  PDFs. 

From  the  physical  point  of  view,  the  model  predictability  is  also  clearly  quantified  by  r-  CDF, 
which  is  the  probability  that  a  model  is  able  to  predict  with  accuracy  higher  than  a  ,  at  times  larger 
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than  t  (see,  also  Section  1).  Accordingly  to  Bury  (1999)  CDF  for  Weibull  distribution  (18)  can  be 
written  as 


P(t)  =  exp  - 


The  parameter  /?  is  a  measure  of  how  quickly  the  r  -  CDF  decays  with  t.  The  super-exponential 


(y#  >1),  exponential  (/?  =  !)  and  sub-exponential  (/?<!)  decay  regimes  correspond  to  small, 


intennediate  and  high  probabilities  to  obtain  the  enhanced  predictability  in  an  individual  forecast. 

The  following  method  was  used  to  estimate  distribution  parameters  of  (18)  from  ensemble  samples. 
The  PWMs  can  be  estimated  in  two  ways:  from  Eq.  (Dl),  using  analytical  representation  (19),  and 
directly  from  ensemble  samples  through  Eq.  (D2).  The  distribution  parameters  are  obtained  through  a 
misfit  between  these  estimations.  For  details  see  Appendix  D. 

Typical  r-CDFs  computed  for  both  WNLPs  and  RNLPs  are  shown  in  Fig.  11a  and  Fig.  lib, 
respectively.  A  more  tailed  CDF  corresponds  to  the  RNLPs  (compare  Fig.  11a  and  Fig. lib).  That 
indicates  a  higher  probability  of  long  individual  forecasting  for  spatially  correlated  initial  errors  than  for 
WNLPs.  If  parameters  for  distribution  (18)  are  known,  there  is  a  chance  to  predict  an  asymptotic 
behavior  of  r  -  CDF  tail  with  t,  and  understand  contributions  of  rare  long  forecasts  to  model 
predictability. 

The  last  statement  is  illustrated  through  a  simple  example.  Let  us  estimate  a  contribution  of  rare 

— 9 

long  forecasts  to  the  growth  of  r  -  mean  with  R  .  The  mean  IPT  calculated  for  WNLPs  (/  =  0.1  ) 


_2 

added  to  initial  conditions,  was  equal  to  about  30.41  days,  when  the  tolerance  level  was  s  =  0.2  .  For 
initial  RNLPs  with  the  same  intensity  but  the  correlation  radius  R_  =  125  km,  the  Weibull  distribution 
parameters  are  estimated  as 

y  «  37.5  days ,  7  «  6.2  days  ,  and  /?  «  2. 1 .  (20) 


For  Weibull  distribution  (18)  the  r  -  mean  is  calculated  as 
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(z)  =  y  +  rir(\  +  \l  p). 


(21) 


where  Y  is  the  gamma  function. 

Substituting  parameters  (20)  into  Eq.  (21)  yields, 

(t)  ~  44.2  days  .  (22) 

Therefore,  the  mean  IPT  has  grown  on  13.8  day.  The  tail  lengthening  stimulates  approximately 
50%  (rj  x  6.2  days)  of  this  growth.  The  other  half  of  the  growth  estimated  as  y  -  (r)  =  7. 1  days  is  caused 
by  variations  of  the  location  parameter. 

9.  Transition  From  Linear  to  Nonlinear  Predictability  Regime 

The  above  results  showed  the  existence  of  two  predictability  regimes.  From  the  practical  point  of 
view  it  is  important  to  understand  how  quickly  the  linear  predictability  decays  and  what  is  a  threshold  8 
on  amplitude  of  initial  perturbations,  above  which  the  nonlinear  predictability  regime  dominates.  There 
are  two  approaches  for  determining  the  duration  of  the  linear  predictability  regime.  One  is  to  compare 
the  evolution  of  a  perturbation  under  the  full  nonlinear  model  with  its  evolution  under  the  tangent  linear 
model  in  order  to  quantify  the  PE  in  this  model  as  a  function  of  time  (for  example  see  Vukicevic,  1991). 

Another  approach  is  to  develop  a  criterion.  For  example  Gilmour  et  al.  (2002)  suggested 
estimating  duration  of  the  linear  predictability  regime  by  monitoring  the  evolution  of  twin  perturbations 
under  the  full  nonlinear  model.  Our  idea  is  to  demonstrate  that  the  decay  of  the  linear  predictability 
regime  correlates  with  the  changes  in  behavior  of  t  -  variance. 

Let  us  introduce  the  non-dimensional  parameter  /u  =  s2  / / 2  that  is  a  measure  for  degree  of 

smallness  of  initial  perturbations,  and  estimate  it  for  a  number  of  combinations  of  s  "  and  /  .  The 

_ •) 

linear  predictability  regime  occurs  for  2.5x10  <jU<  1.0  (Fig.  12  a,  b)  when  r-mean  and 
r  -  variance  evolve  along  the  following  logarithmic  laws: 

(r)  ~  ln^2  / 72),  (23) 
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(24) 


(£r2)~ln (£2 /I2). 

Eq.  (23)  coincides  with  the  well-known  law  of  the  exponential  growth  of  infinitesimal  initial 
perturbations  on  a  chaotic  attractor  (Lorenz,  1996).  Logarithmic  law  (24)  was  theoretically  predicted  by 
Chu  et  al.  (2002). 

The  departure  from  the  linear  predictability  is  observed  for  1.0  <  ju  <  16.7  (Lig.  12c, d). 
Although,  here,  the  r  -  mean  still  follows  the  [quasi]  logarithmic  law  (23)  with  the  growth  (reduction) 

of  ju(I  ),  and  deviations  from  this  law  do  not  seem  to  be  essential  for  small  /  and  large  s  ,  the 
z  -  variance  does  not  demonstrate  a  behavior  that  Eq.  (24)  predicts.  It  is  practically  constant  (about  9 

days  )  for  /  varied  between  0.5x10  and  0.2  (Lig  12d).  That  shows  no  correlations  between  the 
z  -  mean  (ensemble  mean)  and  z  -  variance  (ensemble  spread)  when  the  amplitude  of  initial 
perturbations  exceeds  the  threshold  //  «  1.0 . 

The  intennediate  domain  between  the  linear  and  nonlinear  predictability  regimes,  when  there  are 
no  correlations  between  ensemble  mean  and  spread,  is  also  characterized  by  small  changes  in  scale  and 
shape  parameters  (  r/,/3)  estimated  as 

7«9.6±0.1,  P « 3.4±0.1.  (25) 

Here,  z  -  statistics  is  still  close  to  Gaussian  because  the  skewness  and  the  kurtosis  are  close  to  0 

(|SK|<0.5)  and  3  (|KU-3.0|<  1 .5  x  10_1),  respectively. 

The  lack  of  statistically  significant  correlations  between  the  forecast  skill  and  ensemble  spread  is 
sometimes  observed  in  ensemble  modeling  for  the  atmosphere  and  ocean.  Lor  example,  Moore  (1999) 
performed  a  series  of  experiments,  in  which  he  applied  different  methods  of  ensemble  generation  to  a 
quasi-geostrophic  model  of  the  Gulf  Stream.  He  obtained  no  statistically  significant  relationships 
between  forecast  skill  and  ensemble  spread  in  a  number  of  cases  and  explained  them  only  by  poor 
statistics  in  the  forecast  experiments.  Our  computations  show  that  the  same  effects  can  appear  when  the 
initial  perturbations  have  quite  large  amplitudes. 
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10.  Predictability  Horizon 


Asymptotic  behavior  of  the  v  -  CDF  with  time  determines  the  predictability  horizon  for  the 
hydrodynamic  model  used  here.  Model  predictability  horizon  is  defined  as  the  maximum  prediction  time 
{Thor)  f°r  the  given  model  and  statistics  of  initial  perturbations  (Kravtsov  1993)  and  can  be  calculated 
by 


thor  =7  +  11 


- * 

In  P 


Up 


(26) 


where  P  is  the  probability  that  T/Wr  will  be  achieved  in  an  individual  forecast. 

Let  us  estimate  the  model  predictability  horizon  for  the  example  discussed  in  Section  8. 
Substituting  parameters  (20)  into  Eq.  (26)  yields, 

Thor  *50.21  ,  52.79  and  55.03  days  (27) 


for  P*  =  0.01 , 0.001  and  0.0001,  respectively.  Thus,  for  the  RNLPs  with  the  fixed  values  of  ffy  and 

—  2  _2 

/  ,  and  s  =0.20,  the  model  predictability  horizon  is  limited  to  52-55  days  and  any  individual  forecast 
longer  than  55  days,  is  improbable. 

11.  Conclusions 


We  have  used  a  simple  shallow-water  model  to  understand  predictability  of  perfect  ocean  models. 
This  model  is  a  highly  idealized  representation  of  some  aspects  of  the  ocean  dynamics  and  naturally, 
cannot  simulate  re-distributions  of  PE  between  barotropic  and  baroclinic  dynamics,  and  interactions 
among  large  and  small  flow  scales  reproduced  by  high-resolution  ocean  models.  However,  due  to  the 
small  degree  of  freedom  of  the  model,  distribution  function  for  a  prediction  scale  and  its  high-order 
moments  were  computed  for  a  large  number  of  ensemble  realizations  (up  to  1 05).  That  allows  us  to 
qualify  PE  statistics  in  an  accurate  manner. 

Similar  analysis  is  difficult  to  realize  in  full-scale  numerical  forecast  ocean  models  due  to  limited 
computer  resources.  The  full-scale  models  produce  small  sample  forecast  ensembles  and  therefore, 
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cannot  resolve  the  full  complexity  of  PE  statistics.  The  idealized  model  reveals  some  trends  in  behavior 
of  PE  and  produces  useful  knowledge  for  extracting  them  from  small  forecast  ensemble  samples. 

The  following  trends  obtained  in  the  present  study  need  to  be  examined  by  baroclinic  high-resolution 
ocean  models  with  realistic  bottom  topography  and  forcing. 

(1)  In  the  limit  of  zero  horizontal  viscosity  the  model  demonstrates  the  initial  decay  of  spatially 
uncorrelated  and  correlated  perturbations  due  to  nonlinear  bottom  friction.  Since  the  bottom  friction 
plays  a  significant  role  in  the  energy  balance  of  ocean  coastal  currents  (Wunsch  and  Ferrari,  2004),  the 
observed  decay  seems  to  be  an  important  process  that  may  enhance  model  predictability. 

Initial  error  decay  was  noted  in  connection  with  numerical  experiments  in  preparation  for  what 
became  eventually  known  as  the  Global  Weather  Experiment  (Williamson  and  Kasahara,  1971),  and 
more  recently,  in  the  oceanographic  and  meteorological  literature,  by  Brasseur  et  al.  (1996),  Wirth  and 
Ghil  (2000),  Vannitsem  and  Nicolis  (1997),  Snyder  et  al.  (2003)  among  others.  The  initial  error  decay 
found  in  these  models  was  controlled  by  dissipative  processes  parameterized  through  a  hyper-diffusion 
or  usual  diffusion  operator  with  a  horizontal  viscosity.  Our  study  reveals  another  mechanism  of  initial 
error  decay:  PE  decays  due  to  viscosity  damping  by  nonlinear  bottom  friction  and  this  process  is  weakly 
dependent  on  horizontal  viscosity. 

(2)  Statistics  of  the  finite-amplitude  prediction  error  was  found  to  be  close  to  extremum  statistics  - 
Weibullian.  No  general  theory  is  available  that  demonstrates  universality  of  this  result  for  any 
forecasting  model.  However,  there  is  a  number  of  reasons  to  examine  Weibull  statistics  in  large  ocean 
and  atmospheric  models.  First,  extremum  statistics  are  often  observed  to  arise  in  multi-dimensional 
systems,  exhibiting  correlation  over  a  broad  range  of  scales,  leading  to  emergent  phenomenology,  such 
as  self-similarity  and  in  some  case  fractional  dimension  (Boffetta  et  al.,  2002).  Second,  Weibull  statistics 
seems  to  be  a  good  mathematical  tool  for  the  parametrical  estimate  of  PE  distributions  in  small  forecast 
ensembles  and  from  limited  observation  samples.  Preliminary  computations  provided  by  us  support  this 
conclusion  (Ivanov  et  al.,  2007).  Third,  if  the  divergence  of  a  predicted  flow  in  phase  space  is  constant, 
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the  cumulative  distribution  function  P(s,t)  (defined  by  Eq.  5a)  can  be  analytically  estimated  (Ivanov  et 
al.,  2007).  Here,  we  only  give  the  final  results  without  calculations: 

P(e,t)  =  exp(-c0  t/(r)),  (28) 

where  k  is  a  constant  depending  on  a  forecast  model. 

CDF  (28)  corresponds  to  the  exponential  probability  distribution  function,  which  results  from 
Weibullian  distribution  (18)  for  y  =  0 ,  /?  =  1 ,  and  rj  =  c0  /  (r) .  Eq.  (28)  assesses  predictability  for  many 

forecast  models  including  Lorenz  63  model  (Lorenz,  1963). 

(3)  Our  computations  demonstrated  that  the  transition  from  the  linear  to  nonlinear  predictability 
regime  may  be  detected  by  high-order  t  -  moment  behavior.  We  found  the  predictability  regime  where 
the  mean  prediction  time  scale  grew  along  [quasi] -exponential  law,  which  accompanies  the  linear 
predictability,  but  a  behavior  of  variance  for  this  scale  was  abnonnal,  so  there  are  no  statistically 
significant  correlations  between  the  forecast  skill  and  ensemble  spread.  This  result  explains  the  lack  of 
correlations  between  the  forecast  skill  and  ensemble  spread  found  in  a  number  of  the  atmospheric  and 
oceanic  models,  and  may  be  used  to  develop  criteria  for  detection  of  the  transition  from  linear  to 
nonlinear  predictability  regime. 

Although  the  present  study  has  analyzed  only  model  predictability  near  the  equilibrium  ocean 
state,  non-monotonous  behavior  of  ensemble  spread  seems  to  be  a  common  feature  of  forecast  models. 
A  possible  explanation  of  this  effect  is  that  model  predictability  strongly  depends  on  the  spatial 
correlation  scale  of  initial  error  as  it  was  firstly  demonstrated  by  Lorenz  1965,  and  on  variations  of  this 
scale  when  the  PE  evaluates. 

There  are  in  general  different  scenarios  in  behavior  of  the  ensemble  spread  due  to  variations  of  the 
spatial  correlation  scale  of  PE.  For  example,  Lorenzo  et  al.  (2003)  demonstrated  decreasing  ensemble 
spread  of  an  ensemble  system  consisting  of  chaotic  Lorenz  cells  diffusively  coupled,  as  the  correlation 
scale  of  PE  grows.  If  the  growing  perturbations  rapidly  adopt  the  horizontal  scales  comparable  to  that  of 
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the  reference  flow,  the  ensemble  spread  reaches  an  asymptotic  value  and  then  weakly  changes.  This  case 
is  studied  in  Section  9. 

However  independently  on  the  realized  scenario,  the  correlation  scale  of  prediction  error  should 
change  when  the  error  leaves  the  tangent  space  (Schertzer  and  Lovejoy,  2004).  That  can,  in  general, 
result  into  variations  of  high-order  statistics  of  IPT  and  in  the  non-monotonous  behavior  of  ensemble 
spread,  examples  of  which  can  be  found  in  Nicolis  (1992),  Cohn  (1993)  among  others.  Obviously,  a 
more  systematic  analytical  description  of  the  observed  effects  is  required.  We  will  address  this  point 
elsewhere. 
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Appendix  A 


Spectral  Decomposition 

Following  Morse  and  Feshbach  (1953),  the  geostrophic  stream  function  is  decomposed  into 

M 

^  ~  ^hom  ^ harm  >  ^hom  (x>  0  —  (A1) 

m= 1 

where  ^harm  is  the  hannonic  function  calculated  with  the  open  boundary  conditions  written  by 

x2 

^ harm  \  A'  =  ~  J  Uh  (*,  t)dx  ,  and  V\nm  \A=0.  (A2) 

la 

The  basis  functions  {  xVm  }  are  the  eigenfunctions  of  the  horizontal  Laplacian  operator  A  and 
computed  by 

^  W m  ~  m  ’  W m\  AuA'  ~  (A3) 

where  Am  are  its  eigenvalues. 

Spectral  decomposition  (A1)-(A3)  is  also  applicable  for  non-rectangular  domains  and  can  be 
generalized  for  2D  compressible  and  3D  incompressible  flows  (Chu  et  al.  ,  2003). 

The  optimal  truncation  number  M  t  in  (Al)  depends  on  complexity  of  a  flow  structure.  Direct 

computations  have  shown  that  to  represent  the  reference  flow  with  accuracy  better  than  0.1%,  the 
truncation  number  M  t  should  be  taken  about  10". 

This  number  can  be  theoretically  confirmed  in  the  following  way.  The  used  hydrodynamic  model 
may  resolve  ocean  flows  with  -100-200  km  spatial  scales  ( L  )  because  of  numerical  grid  with  50  km 
cells.  Following  Mikhlin  (1964)  the  spatial  scale  of  the  highest-order  mode  y/m  can  be  defined  as 

lM  ~  L\  (A  ^m)V2,  (A4) 

where  L\  is  the  spatial  characteristic  scale  of  the  largest  mode  m  =  1 . 
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8  7  7 

For  the  rectangular  domain  A  we  can  suppose  L  «  L ^  ,  L\  «  10  km  and  Z  « 10  -2x10  km. 


Substituting  L\  and  L M  into  Eq.  (A4),  yields 


-25-100,  (A5) 


Therefore,  from  80  to  120  modes  should  approximate  the  reference  and  perturbed  solutions  with  a 
reasonable  accuracy.  We  used  M opt  =  100  . 
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Appendix  B 


Dissipation  of  Prediction  Error 

Following  Blumberg  and  Mellor  (1987),  the  bottom  friction  of  both  velocity  components  is 
parameterized  by 

Fl(ul,u2)  =  ccEul,  (Bl) 

F2(ux,u2)  =  ccEu  2,  (B2) 

where  a  is  the  drag  coefficient,  E  =  ■s]ul  +u2  .  Let  the  circulation  in  a  basin  be  decomposed  as  the 
reference  circulation  (ux , u2 )  and  perturbations  (u[,u'2)\ 

ul=ul+u[,  u2=u2+u'2, 

such  that 

|h"j  |  »  | u[  |  and  | u2 1 »  \u'2 1 . 

For  small  perturbations,  linearization  of  (Bl)  and  (B2)  leads  to 

F)  m  Fl(ul,u2)  +  aE~112  [(2mj2  +w22jw1'  +ulu2u'2\  ,  (B3) 

F2  ~  F2  (ul ,u2)  +  aE~l/ 2  [(m | 2  +  2u22  ]b2  +  ux  u2u[  ],  (B4) 

Therefore,  the  dissipative  terms  in  the  prediction  error  energy  balance  can  be  written  by 


Fxii[ )*ccE  1/2  Y\ ({u[ )2 )  +  y2 {u[ u2) 


F2u2)  «  aE~112 


f3((w2)2)  +  r2(wi'“2). 


(B5) 


(B6) 


where  yx  =2u2  +u22  ,y2  =ulu2  and  y3  =u2  +2  u2  . 
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Appendix  C 


Probability  Density  Function  of  Initial  Perturbations 

Let  the  two-point  correlation  function  of  the  geostrophic  stream  function  be  projected  onto  the 
phase  space  with  the  basis  [  ]  , 


M  M 


)^(x2  )> =2^2^  \dA'dA™ )  ^ (X2 } 


(Cl) 


l=\  m= 1 


Integrating  (Cl)  over  the  spatial  variables  x1;x2  in  the  computation  area  yields 

^d\d2x2{SW(xl)SF(x2))  =  (SA?BlmSA°m),  (C2) 


with  jjVx1</2x2'F/  (x,  )H>m  (x2 )  =  Bhn . 


(C3) 


Using  (10),  the  correlation  function  in  the  left-hand  side  of  (C2)  is  represented,  following  Panchev 
(1971)  as 

\2  \ 


(8'F(tlx)S'F(tl2))  =  1/2 Rif-  exp 


(xi  -XiY 

v  2R^  j 


R±  «  L  ~  L, . 


(C4) 


Introducing  new  variables:  r  =  xx  -  x2  and  r,  =  x ,  +  x 2 ,  and  integrating  (C2)  over  r  and  r,  yield 


SA !,XX)  =  2k  LiL22R,J2ern-jF-)erf( -jSL-) . 


(C5) 


Eq.  (C5)  defines  an  M-dimension  hyper  ellipsoidal  surface  Qq  in  the  model  phase  space.  The 
ellipsoid  semi-axes  are  directed  along  the  eigenvectors  of  the  matrix  Blm .  Therefore,  the  maximum 


deviation  for  the  vector  of  initial  perturbations  ((8A^ax)2  }  is  estimated  by 


1/2 


((*C)2)''2  =  (2k)V2l1l2r±2I 


y/2R ,  '  y[2R 


(C6) 


where  2niax  is  the  maximum  eigenvalue  of  the  matrix  Blm  . 
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In  the  context  of  the  predictability  problem,  the  PDF  of  a  RNLP  is  Gaussian  with  the  mean 


SA®n  )  =  0  and  the  covariance  matrix 


<rfm  =  BlmLfL?R;4r2 


’lm^\ 


erf( 


L, 


U 


-)~l  erf (~j=^  )_1  . 
V2 R ,  V2 R , 


(C7) 
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Appendix  D 


Probability  Weighted  Moments 

The  probability- weighted  moments  of  a  random  variable  are  defined  by  Greenwood  et  al.  (1979), 
1  1 


ak 


X (u)(\  -  u)k  du  ,  and  Afc=J 


X(u)u  du  , 


(Dl) 


0 


0 


where  k=0,l,...,K;  X(u)  is  the  quantile  function  (  i.e.,  the  inverse  of  cumulative  distribution  function). 
Following  Hosking  and  Wallis  (1997)  (CX/(,  /4  )  are  calculated  from  an  ordered  random  sample  rn  of 

size  N  as  unbiased  estimates  of  (ak  ,bk): 


1  n 

a  i,  = — C 
N  1 


'N- 1^1  — 


z 


c 


N  -  ti\ 


V  k  )  ,  V  k  J 

n= 1 


Tn  >bk 


f N-\ '~X  — 
k 


X 

n= 1 


c 


n  —  1 


rn ,  (D2) 


where  C 


are  the  binomial  coefficients. 


The  following  procedure  is  used  to  estimate  the  parameters  of  a  Weibull  distribution  function. 

(a)  The  probability  weighted  moments  a/(  or  /4  are  computed  from  modeled  samples  accordingly  to 

(D2).  (b)  With  the  given  moments  (Ct/c ,  J3k)  ,  the  Weibull  distribution  parameters  are  identified 
from  the  condition: 


KLS  (j 7 ,  p,  y)  =  KL(Fl  \F2  )  +  KL(F2  |F, )  ->  min  ,  (D3) 

where  KLS  is  the  symmetrical  Kullback-Leibler  distance,  F\  and  F2  are  distribution  functions 
computed  from  Monte-Carlo  samples  and  identified  by  the  probability  weighted  moments  ,  respectively. 
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T-t  Tz  t  (days) 

Fig.  1  The  IPT  and  e-folding  time  for  an  oscillating  prediction  error,  (a)  IPT  ( r  )  computed  in 
an  individual  forecast.  Shaded  zones  show  returns  of  model  predictability,  (b)  an 
ensemble  of  J  (dashed  curves),  the  ensemble  averaged  J  (solid  curve)  and  the  e-folding 
time  ( ze  ),  and  (c)  the  ensemble  of  J  (dashed  curves)  and  an  appropriate  ensemble  of 

IPT( 
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Fig.  2  Two  areas  for  POM  integration:  Domain-A  with  the  rigid  A  (the  boundary  segment  between  lb  and  la  in  the  counter¬ 
clockwise  direction)  and  open  A'  boundaries  (top  panel),  and  Domain-B  with  the  rigid  boundary  (bottom  panel). 
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Fig.  4  Sensitivity  of  z  -  statistics  to  the  ensemble  size  (N):  (a)  the  Kullback-Leibler  distance, 
(b)  the  z  -  mean  ;  (c)  the  z  -  most  probable  ;  (d)  the  r  -  variance;  (e)  the  z  - 
skewness  and  (g)  the  z  -  kurtosis.  Initial  perturbations  are  red  noise  with  the 

—2  —2 

correlation  radius  =  1 12  km  and  the  intensity  /  =0.1,  s  =  0.2  .  Arrows  indicate 

^=103. 


Fig.  5  Temporal  evolution  of  the  root  mean  square  error  (j'j  for  various  initial  perturbations: 

(a)  WNLPs  with  different  noise  intensities  (I2):  0.05  ,  0.01,  0.005,  and  0.001 
(denoted  by  ‘1’,  ‘2’,  ‘3’,  ‘4’),  and  (b)  RNLPs  with  the  correlation  radius  ( R j_)  of  70 

km  and  different  noise  intensities  (I2):  0.02  ,  0.01,  0.003  and  0.001  (denoted  by  ‘1’, 
‘2’,  ‘3’,  ‘4’). 
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Fig.  6  The  growth  rate  Q  for  (a)  NWLPs  and  (b)  RNWPs  with  the  correlation  radius  (  Rj_  )  of  70  km. 

— 9  —9 

The  solid  and  dashed  curves  correspond  to  /  =1.0x10  and  /  =0.5.  Black  arrows  indicate 
boundaries  between  the  linear  and  nonlinear  predictability  regimes  for  different  characteristics  of 
initial  perturbations. 
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Fig. 8  The  spectral  entropy  S  of  (a)  PE  and  (b)  the  reference  solution.  Dashed  and  solid  lines  are 
computations  for  the  first  20  ( m0  =1,M  =  20)  modes  and  last  50  ( mQ  =50,M  =  100)  modes, 

respectively.  The  initial  error  is  the  WNLP  with  I  =0.1.  f  =0.5. 


Fig.  9.  t  -  PDFs  for  (a)  WNLPs  with  / 2  =0.1  ,  and  a2  =  0.5  ,  and  (b)  RNLPs  with  I2  =0.1 
and  Rj_  =  125  km ,  and  s  =  0.2  .  Skewness  and  kurtosis  are  0. 15  and  3.09  in  case  (a), 
and  0.77  and  3.95  in  case  (b).  Dashed  lines  indicate  mirror  reflections  of  the  left  hand 
side  tails  of  r  -  PDFs. 
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Fig.  10.  Dependence  of  z  -  statistics  on  the  correlation  radius  Rj_  for  RNLPs  with  the  noise 

—  D  _D 

intensity  /  =0.1  ,  and  s  =0.2.  (a)  r-mean,  (b)  r  -  variance,  (c)  z  -  skewness,  (d) 
z  -  kurtosis. 
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Fig.  11.  t  -  CDF  computed  for  (a)  WNLPs  with  7  2  =0.1  and  (b)  RNLPs  with  7  2  =0.1  and 
R j_  =  70  km  .In  both  these  cases  a1  =0.1.  Here,  the  circles  and  dashed  lines  represent 
t  -  CDFs  computed  directly  from  the  Monte-Carlo  samples  and  by  the  probable 
weighted  moment  technique,  respectively.  The  location  parameter  y  equals  to  20.4  days 
in  case  (a)  and  38.3  days  in  case  (b). 
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Fig.  12.  Dependence  of  r  -  statistics  on  the  noise  intensity  for  WNLPs:  (a)  and  (b)  for  s  =  0.01  ; 

—2 

(c)  and  (d)  for  e  =0.5  .  Logarithmic  law  (23)  is  indicated  by  solid  line  in  (a)  and  (c). 

Dashed  line  in  (b)  corresponds  to  logarithmic  law  (24). 
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